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ABSTRACT 


The  collapse  of  a  homogeneous  fluid  mass  immersed  in  a  stably 
stratified  fluid  is  studied  numerically.  A  finite  difference  formu¬ 
lation  of  the  Navier-Stokes  equations  in  the  primitive  variables  is 
solved  in  a  large  box  several  times  the  size  of  the  mixed  region. 

The  formulation  conserves  total  energy  in  the  box  in  the  special  case 
where  the  viscosity  is  zero.  The  shape  of  the  homogeneous  region  and 
its  energy  content  are  followed  in  detail .  Confirming  a  previous 
speculation  made  from  a  crude  analytical  theory,  most  of  the  energy 
in  the  homogeneous  fluid  mass  is  shown  to  be  transferred  to  the 
exterior  fluid  in  one  Brunt-Vaisala  period.  The  predictions  agree 
with  available  analytical  models  in  initial  and  intermediate  stages 
and  with  a  previous  tank  experiment  in  the  intermediate  and  late 
stages  of  colxapse. 
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Introduction 

The  flow  phenomena  accompanying  the  collapse  of  a  homogeneous 
fluid  inass  immersed  in  a  stratified  fluid  have  been  the  center  of 
some  attention  over  the  past  ten  years.  The  fluid  mechanical  model 
has  various  applications  in  geophysics  and  in  ';ngineering  but  the 
flow  phenomena  are  complex  enough  to  have  yielded  only  marginally  to 
analysis . 

The  model  is  taken  to  be  a  cylmdr ically  shaped  homogeneous 
fluid  mass  immersed  in  a  linearly  stratified  fluid.  The  cylindrical 
region  is  initially  circular  in  shape  and  the  density  of  the  mixed 
fluid  within  is  equal  to  the  density  of  the  exterior  fluid  at  the 
axis  of  the  cylinder.  Since  the  surrounding  fluid  is  assumed  to  be 
stably  stratified,  the  fluid  in  the  upper  part  of  the  mixed  region  is 
heavier,  and  that  in  the  lower  part  lighter,  than  its  surroundings. 
Thus,  when  released  from  rest,  the  region  expands  horizontally  while 
collapsing  vertically  to  the  level  of  equilibrium  of  the  homogeneous 
fluid.  The  motion  of  the  mixed  region  occurs  at  the  expense  of  the 
potential  energy  originally  stored  in  the  initial  configuration. 

Tht.  distortion  from  the  circular  shape  causes  motions  in  the  exterior 
fluid  that  eventually  radiate  the  energy  away  in  the  form  of  internal 
gravity  waves. 

Previous  studies  aimed  directly  toward  the  collapse  process  have 
been  mostly  experimental  (cf.  Schooley  and  Stewart  L 1 J »  Schooley  .2^, 
Wu  5]  )  .  There  have  been  several  numerical  experiments  Vessel  „ 
Padmanabhan  et  al  [5],  Vasiliev  et  al  '6])  but  only  one  relevant 
analytical  effort  (Bell  and  Dugan  7]).  There  have  been  analytical 
contributions  to  certain  aspects  of  the  collapse  process  by  Hartman 
and  Lewis  [B 1  and  Schooley  and  Hughes  Jl  but  these  will  be  excluded 
from  discussion  because  of  severe  limitations  of  the  linearize! 
equations  of  motion  in  the  fully  mixed  case. 

Numerical  experiments  aimed  at  clarifying  the  physical  processes 
should  be  rrost  enlightening  but  they  have  not  been  entirely  success-  ^ 
f ul .  The  physical  model  solved  numerically  by  Padmanabhan  et  al 
is  an  idealized  one  in  which  fluid  motions  exterior  to  the  homogeneous 
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region  are  ignored.  Thus,  no  energy  can  be  transferred  to  the 
exterior  fluid  and  the  model  could  be  an  accurate  one  only  for  the 
initial  stage  of  collapse.  In  the  initial  stage,  though,  those 
predictions  are  at  variance  with  an  analytical  soluLion  of  the  same 
model  (cf.  Bell  and  Dugan  [7j )  and  it  appears  that  the  solution  did 
not  conserve  energy.  Wessel  [1*]  solved  the  Navier-Stokes  equations 
numerically  in  such  a  way  that  exterior  fluid  motions  were  included 
in  the  computation.  Gross  corroboration  with  the  experiment  of 
Wu  [,5l  was  obtained  but  the  calculation  was  not  accurate  enough  to 
attempt  detailed  predictions  of  energy  balances.  Finally,  Vasiliev 
et  al  6]  have  solved  numerically  a  particular  physical  model  that 
retains  some  aspects  of  the  growth  of  a  mixed  region  due  to  turbulence 
as  well  as  the  resulting  collapse.  However,  no  comparisons  with  other 
data  or  conclusions  were  made. 

In  summing  up  the  previous  contributions,  there  remain  several 
outstanding  questions  about  the  mechanics  of  the  collapse  process. 

First,  the  results  for  the  initial  stage  of  collapse  are  contradic¬ 
tory.  Bell  and  Dugan  "f\  discuss  several  analytical  models  for  the 
initial  stage  that  would  appear  to  be  mechanically  convenient.  How¬ 
ever,  in  those  models  the  homogeneous  region  is  assumed  to  conserve 
its  energy  and  this  is  a  questionable  hypothesis.  The  experiments 
of  Wu  i]  are  not  helpful  in  answering  this  question  because  the 
experimental  method  masked  the  initial  stage,  and  the  numerical 
predictions  of  Wessel  [U]  and  Vasiliev  et  al  [6] are  not  accurate 
enough.  Second,  there  has  been  no  analytical  or  numerical  method  that 
Is  adequate  to  predict  the  phenomena  in  the  late  stage  of  collapse. 

Last  of  all,  the  energetics  of  the  motions  and  the  resulting  physical 
implications  have  not  been  explored. 

This  paper  is  written  to  answer  these  questions.  The  Navier- 
Stokes  equations  in  the  Boussinesq  approximation  are  solved  numerically 
in  a  large  box.  The  fluid  is  assumed  to  have  constantly  increasing 
density  with  depth  except  for  a  circular  region  that  is  initially 
constrained  to  be  homogeneous.  The  numerical  method  follows  along 
lines  laid  out  by  Williams  l'J  and  Piacsek  and  Williams  1L  and  it 
conserves  the  total  energy  in  the  box  remarkably  well.  The  source 
of  energy  for  the  fluid  motions  is  the  potential  energy  initially 
stored  in  the  homogeneous  region.  The  rate  of  energy  transfer  from 
the  homogeneous  region  to  the  surrounding  fluid  is  followed  in  detail 
and  the  result  substantiates  the  speculation  of  Bell  and  Dugan  7] 
that  the  mixed  region  does  not  conserve  its  energy  even  in  the  early 
stage.  Also,  confirming  another  speculation  of  Bell  and  Dugan  Tj , 
the  viscous  stresses  are  shown  to  affect  the  solution  only  at  a  late 
stage  of  collapse . 

Numerical  Analysis 

The  Navier-Stokes  equations  are  solved  by  an  adaptation  of  a 
method  proposed  by  Williams  CIO]  and  Piacsek  and  Williams  Ill]  .  The 
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equations  of  motion  are 


(1) 


where  (u, w)  are  the  (  x,  z)  components  of  velocity,  p  is  the  pressure, 
T  is  the  temperature,  0  is  the  kir emetic  .iscosity,  and  ^  is  the 
thermometric  diffusivity.  The  operator  is  the  convective 

der i  /ative  and  the  density  is  related  to  the  temper¬ 
ature  by  the  relation  {  =  1  -  *(.T“T0^}  where  <  is  the  thermal 

expansion  coefficient.  The  first  two  equations  represent  the  con¬ 
servation  of  momentum  within  the  Roussinesq  approximation  wherein 
density  changes  are  assumed  negligible  ir,  the  acceleration  terms. 

The  third  equation  represents  the  condition  of  incompressibility  and 
the  last  one  the  conservation  of  energy.  These  equations  are  to  be 
solved  subject  to  approximate  boundary  and  initial  conditions.  Tile 
boundary  conditions  are  taken  to  be  no-stress  conditions  with  the 
temperature  (density)  fixed  at  the  top  and  bottom  of  a  rectangular 
box.  In  detail,  the  boundary  conditions  are 
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where  T  is  fixed  on  z  =  o  and  d,  and  b  and  d  are  the  horizontal 
and  vertical  sizes  of  the  box,  respectively.  The  initial  condition 
is  specified  as  quiescent  with  horizontal,  equally  spaced  isopycnals 
except  for  a  circular  region  in  which  the  fluid  is  totally  mixed. 
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A  Poisson  equation  for  the  pressure  can  be  derived  from  equations 
(l)  so  that 

£  +  *•§  V 

vhere  is  the  velocity  divergence  and  g  represents  the  nonlinear  and 
viscous  terms  in  the  Navier-Stokes  equations.  The  integral  form  of 
this  equation  is 

*o 

so  that 

"O  (4) 

lias  to  be  satisfied  around  the  boundary.  Since  the  no-stress  boundary 
conditions  are  to  be  imposed,  the  boundary  condition  on  equation  (3) 
is  of  Neumann  type  and  it  is  important  that  the  finite-difference  form 
of  equation  (3)  satisfy  the  integral  (4)  on  the  boundaries. 

In  order  to  establish  the  finite-difference  scheme  for  solving 
the  above  equations,  the  difference  operator 

U -4  <*-¥>] 

and  the  sum  operator 

are  defined  in  the  notation  of  Grammeltvedt  f  12]  .  These  operators 
form  a  linear  commutative  and  distributive  algebra  for  which  various 
operator  rules  and  identities  can  be  constructed,  such  as 

A  quadratic  ally  conservative  scheme  in  the  sense  of  Piacsek  and 
Williams  (11]  can  be  constructed  for  the  rele.avt  equations  as 

**T  *  =  -ix(.^T  *)  -  S*(«-T  ')  *■  iT( +  l^) 

•.  K.  <.i„T  +  VjO 

) 
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where  variables  are  defined  on  the  staggered  grid  as  snown  in  Figure  1. 


Figure  1 


Computation  grid 

•  horizontal  velocity 
X  ■  vertical  velocity 

^  X  pressure  and  temperature 


The  terms  of  the  form 

force  the  finite-difference  scheme  for  the  nonlinear  terms  to  be. 
quadratically  conservative  regardless  of  whether  the  divergence  is 
zero  or  not.  In  other  words,  the  quadratic  conservation  is  algeberaic 
and  is  independent  of  the  accuracv  of  the  solution  in  the  sense  of 
how  close  the  divergence  is  to  being  zero.  This  does  introduce  an 
error  proportional  to  the  divergence  (.which  is  not  identically  zero 
during  the  computation)  into  the  integrals  of  linear  quantities. 
However,  these  linear  conservation  integrals  are  not  . as  necessary  or 
as  meaningful  a  requirement  for  computational  stability  as  are  the 
quadratic  ones . 

The  nonlinear  terms  are  evaluated  at  time  t,  thus  constituting 
the  "leap-frog"  method  (Richtmyer  [13],  p.  17)  which  has  a  time 
truncation  error  of  0(  t2)  and  a  von  Neumann  condition  of 

[  H  *  t%  # 

The  viscous  terms  are  lagged  in  time  at  time  level  t-*t  and  are  subject 
to  the  stability  criteria  of 

*  <  ^>6, 

and  the  pressure  and  temperature  terms  are  time  centered  (evaluated 
at  time  t  ) . 

Another  instability  that  can  appear  is  aliasing  (Phillips  [1;0 ) . 
This  instability  arises  when  waves  that  are  too  short  to  be  resolved 
by  a  given  set  of  grid  points  are  misrepresented.  Such  instabilities 
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are  controlled  s\x:cessfully  by  using  schemes  like  this  one  that  main¬ 
tain  the  integral  constraints  of  physical  importance  on  the  quadratic 
quantities  u2,  w2,  and  T2.  For  a  three-level  time  scheme  like  this 
one,  computational  modes  are  present  due  to  the  fact  that  a  first  order 
continuum  equation  in  time  has  been  raised  to  a  second  order  difference 
equation  in  time.  The  computational  modes  are  usually  small  and  it 
takes  several  hundred  time  steps  for  serious  deviations  from  the 
correct  solution  to  arise.  These  computational  modes  are  suppressed 
by  periodically  averaging  over  adjacent  time  steps. 


The  Harlow  and  V/eleh  [15]  numerical  strategem  to  drive  the  civergence, 
<s>,  to  zjro  was  applied.  This  is  necessary  because  a  degree  of  round¬ 
off  error  is  inevitable  and  it  leads  to  the  creation  of  an  artificial 
divergence.  The  computational  strategem  works  as  follows:  since  the 
divergence  ii^the  computation  at  a  given  time  step  is  not  exactly 

zero  (  4^o)  hut  tha',  of  time 't+d't  ought  to  be  zero 

is  treated  as  equal  to  -  eh  . 

The  actual  solution  of  equation  (6)  with  the  predictors  (5)  is 
by  partial  Fourier  reduction.  The  boundary  conditions  (2b)  on  the 
pressure  are  made  homogeneous  by  a  simple  transformation  and  the 
pressure  is  then  expanded  in  a  Fourier  cosine  series  in  the  vertical 
direction.  The  expansion  and  decomposition  are  accomplisned  by  means 
of  fast  Fourier  transforms.  The  equation  that  has  to  be  solved  then 
takes  the  form 

^  OK  *  (7) 

where  n  refers  to  the  Fourier  mode  number,  is  the  eigenvalue  for 

the  nth  mode,  and  is  the  transformed  source  function.  This  is  a 

trid iagonal  equation  in  x  and  it  can  be  solved  subject  to  the 
Neumann  boundary  conditions  by  means  of  a  tridiagonal  algorithm 
(Varga  rl61).  For  the  degenerate  case  n=0,  the  matrix  is  no  longer 
liagmally  dominant  and  the  algorithm  breaks  down  so  the  equation  is 
solved  as  a  marching  problem  in  this  case.  The  integral  constraint 
(4)  is  satisfied  by  the  finite-difference  formalism  within  round-off 
so  this  degenerate  case  poses  no  problem. 

The  normal  checks  on  the  computation  for  consistency  and  conver¬ 
gence  were  made  but,  since  the  computation  of  Padmanabhan  et  al  [5] 
apparently  dissipated  energy  in  a  physical  model  where  none  should 
have  been  dissipated,  it  was  felt  especially  important  to  prove  that 
the  present  method  did  conserve  energy  at  least  approximately  in  the 
case  of  zero  viscosity. 
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TIME  INI) 

Figure  2 

Total  energy  in  computation  box 


In  Figure  2,  the  total  energy  in  the  box  is  plotted  versus  time.  The 
time  axis  is  scaled  by  the  Brunt-Vaisala  frequency  N  where 


4t 


and  p  is  the  density  profile  of  the  exterior  fluid.  The  energy  is 
conserved  to  within  5°/0  for  l4  periods  which  is  over  2500  time  steps. 
This  is  remarkable  considering  that  the  initial  condition  data 
includes  a  discontinuity  and  that  in  such  a  long  run  the  mixed  region 
is  highly  elongated  so  that  it  virtually  disappears  into  the  grid 
spac ing . 


Results  and  discussion 


The  box  size  chosen  for  the  computations  is  elongated  so  that 
the  collapsing  mixed  region  can  expand  unrestricted  by  the  sidewall 
boundary.  The  initial  mixed  region  is  taken  +o  be  a  quarter  circle  in 
a  corner  of  the  box  since  the  problem  exhibits  both  horizontal  and 
vertical  symmetry  about  the  mixed  region  centerline.  The  distance  to 
the  top  of  the  box  is  4  radii  and  that  to  the  side  is  l6  radii.-  The 
number  of  grid  points  was  varied  somewhat  to  check  for  convergence 
of  the  solution  but  the  majority  of  the  results  are  for  50X200  grid 
points . 

Figure  5  shows  a  sequence  of  mixed  region  shapes  taken  at 
intervals  of  one-sixth  of  a  Brunt-Vaisala  period.  The  shape  is 
obtained  from  particle  tracers  that  are  massless  points  pushed  along 
by  the  velocity  field.  The  short-time  history  of  the  width  of  the 
mixed  region  is  shown  in  Figure  4  along  with  the  theoretical  pre¬ 
diction  of  Bell  and  Dugan  [7]  for  the  accelerative  stage.  Actually, 
the  numerical  result  is  comprised  of  two  curves  which  bound  the 
region  size  above  and  below;,  the  uncertainty  is  entirely  in  the  plot¬ 
ting  routine  utilized  to  display  the  results.  That  theory  (as  well  as 
similar  ones  of  Padmanabhan  et  al  [51  and  Mei  fiT]  )  essentially 
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Figure  3 

Mixed  region  shape 


Figure  4 

Mixed  region  width 


assumed  that  the  potential  energy  that  was  originally  stored  in  the 
mixed  region  is  converted  to  kinetic  energy  but  no  energy  is  trans¬ 
ferred  to  the  surrounding  fluid.  In  other  words,  that  theory  neglects 
all  motions  in  the  surrounding  fluid.  This  can  only  be  an  arproxima- 
tion  since  all  changes  of  shape  of  the  mixed  region  force  accommo¬ 
dating  motions  in  the  surrounding  fluid.  It  is  plausible  to  extend 
that  theory  slightly  by  making  the  assumption  that,  for  short  times, 
all  motions  inside  the  mixed  region  cause  comparable  motions  outside 
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the  region.  Assuming  that  the  kinetic  energy  of  the  motions  in  the 
two  regions  are  equal,  the  governing  equation  for  the  mi\ed  region 
width  in  the  accelerative  stage  that  is  comparable  to  the  governing 
equation  derived  by  Bell  and  Dugan  [7]  is 

£.(!•* *■  0.  ^ 

The  solution  of  this  equation  is  identical  to  that  of  Bell  , and  Dugan 
[7]  except  that  the  time  axis  is  shifted  by  a  factor  of  2  .  This 

curve  falls  right  on  the  mean  of  the  numerical  results  so  this 
modified  assumption  is  at  least  consistent  with  the  numerical  result. 
This  result  is  interesting  because  of  its  physical  implications. 
Evidently,  the  potential  energy  of  the  mixed  region  initially  is 
lost  to  an  equipartition  of  kinetic  energy  between  the  fluid  inside 
sind  outside  the  mixed  region  and  no  potential  energy  is  stored  in  the 
exterior  fluid  in  this  stage  of  collapse. 

Figure  5  is  another  plot  of  the  width  of  the  mixed  region  versus 
time.  The  dots  are  the  raw  data  from  the  tank  experiments  of  Wu  [3] 


Figure  5 

Mixed  region  width 


and  the  solid  line  is  the  intermediate  stage  theory  of  Bell  and  Dugan 
[7]  .  The  circles  indicate  the  numerical  prediction  of  J he  region 
width  in  the  case  of  zero  viscosity  and  the  triangles  indicate  the 
prediction  in  which  the  molecu- or  viscosity  is  on  une  scale  used  by 
Wu  [3]  .  The  agreement  of  the  numerical  results  with  the  theory  and 
the  experiment  is  quite  good.  This  result  confirms  the  speculation 
of  Bell  and  Dugan  [7]  that  viscous  effects  are  important  only  late 
in  the  collapse.  It  also  sheds  some  doubt  on  the  accuracy  of  the 
numerical  methods  used  by  Paimanabhan  et  al  [5]  since  the  inviscid 
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. . . . 


and  viscous  results  differed  much  earlier  in  time  in  those  predictions. 


Some  physical  insight  is  gained  by  considering  the  energetics 
of  the  motions.  From  the  mixed  regicn  shapes  of  the  form  shown  in 
Figure  3,  along  with  tabulations  of  energy  densities,  it  is  possible 
to  determine  the  amount  of  potential  and  kinetic  energy  in  the  mixed 
region  au  any  given  time.  Figure  7  is  a  plot  o.  these  ei  ergies  versus 


Figure  6 

Energy  balance  of  mixed  region 

"p  ,  potential  energy,  [T] 

K  t  kinetic  energy,  [7] 

5  ,  residual  energy,  [7] 

©  potential  energy,  numerical 
A j  kinetic  energy,  numerical 
0  j  residual  energy,  numerical 

time.  The  energies  are  normalized  by  the  total  potential  energy 
originally  stored  in  the  mixed  region.  These  computations  are  for 
zero  viscosity  (as  shown  above,  viscous  effects  are  not  important  in 
the  early  stage)  so  the  energy  transferred  to  the  exterior  fluid  is 
just  the  total  eneigy  minus  the  sum  of  the  potential  and  kinetic 
energies  in  the  mixed  region.  The  analytical  predictions  of  the 
theories  of  Bell  and  Dugan  [7]  are  also  illustrated  in  Figure  6.  It 
is  remarkable  that  the  results  agree  as  well  as  they  do  because  the 
analytical  results  for  the  energetics  assumed  an  elliptical -shaped 
mixed  region  even  for  times  for  which  the  numerical  results  show  the 
shape  to  be  quite  different  from  an  ellipse.  The  numerical  result 
confirms  the  analytical  prediction  that  most  of  the  energy  has  been 
transferred  to  the  surrounding  fluid  by  the  end  of  one  Brunt-Vaisala 
period  and  it  leads  to  the  conclusion  that  the  energy  transfer 
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mechanism  operates  quickly.  In  turn,  this  leads  to  the  conclusion 
hhat  flow  models  that  assume  energy  conservation  in  the  mixed  region 
(including  those  of  Mei  [17]  ,  Padmanabhan  et  al  [5]  ,  and  the  early 
time  models  of  Bell  and  Drgan  [7]  )  have  very  limited  applicability 
to  the  real  problem  of  mixed  region  collapse. 

A  plot  of  the  potential  and  kinetic  energy  in  the  whole  box 
versus  time  is  also  of  interest  and  is  shown  in  Figure  7.  Although 


Figure  7 

Energy  balance  of  computation  box 


initially  all  of  the  energy  is  potential,  the  collapse  transforms 
roughly  one-quarter  of  it  to  kinetic  energy.  This  prediction  is 
consistent  with  the  fact  that  there  is  an  equipartition  between  kinetic 
and  potential  energy  in  small  amplitude  internal  gravity  waves.  The 
time-averaged  mean  of  the  kinetic  and  potential  energies  is  not  the 
same  but  this  discrepancy  is  attributed  to  the  method  of  computing  the 
potential  energy.  The  zero  level  of  potential  energy  is  assumed  to 
be  that  of  the  linearly  stratified  fluid  outside  the  collapse  region. 
However,  in  a  box  of  finite  size,  the  final  equilibrium  state  is 
distorted  from  the  linear  stratification  by  the  homogenous  fluid  lying 
on  the  mixed  region  axis.  This  distortion  fully  accounts  for  the 
discrepancy  between  the  means  of  kinetic  and  potential  energies  so 
that  equipartition  is  in  fact  achieved. 
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